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Abstract 

A kinetic approach is adopted to describe the exponential growth of a small 
deviation of the initial phase space point, measured by the largest Lyapunov 
exponent, for a dilute system of hard disks, both in equilibrium and in a 
uniform shear flow. We derive a generalized Boltzmann equation for an ex- 
tended one-particle distribution that includes deviations from the reference 
phase space point. The equation is valid for very low densities n, and requires 
an unusual expansion in powers of l/|lnn|. It reproduces and extends re- 
sults from the earlier, more heuristic clock model and may be interpreted as 
describing a front propagating into an unstable state. The asymptotic speed 
of propagation of the front is proportional to the largest Lyapunov exponent 
of the system. Its value may be found by applying the standard front speed 
selection mechanism for pulled fronts to the case at hand. For the equilibrium 
case, an explicit expression for the largest Lyapunov exponent is given and for 
sheared systems we give explicit expressions that may be evaluated numeri- 
cally to obtain the shear rate dependence of the largest Lyapunov exponent. 

KEY WORDS: Lyapunov exponent; shear; Boltzmann equation; pulled 
fronts. 
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I. INTRODUCTION 



Often, familiar models from statistical mechanics exhibit the strong dependence on the 
initial phase space point that we know as chaos . Examples are the Lorentz gas, the hard 
disk and the hard sphere gas 0] . We call a system chaotic if the growth of a deviation from 
a reference trajectory in phase space is exponential, oc exp(A^t), in the limit that the initial 
deviation becomes infinitesimally small. A"*" is called the largest Lyapunov exponent. 

The systems we discuss in this paper, consist of N hard disks with equal mass m and equal 
diameters a, in a two dimensional volume V. The density is n = N/V, and a (dimensionless) 
reduced density is defined as n = na^. The position and velocity of disk / are denoted by 
ri and vi respectively. We define a temperature in equilibrium by NksT = {YliLi 
where the brackets denote an ensemble average. For non-equilibrium stationary states, 
a generalization of this is used for the temperature. A typical velocity vq is defined as 



^^0 = ^JkBT/m. 

Considering a point in phase space {r^, v^} {h = 1 . . . N) and an adjacent point {r^, v^} 
with = rfi + Sr^ and v*^ = Vh + Sv^, one may obtain the largest Lyapunov exponent by 
studying the exponential growth of the deviations Svh and Svh- Other Lyapunov exponents 
exist, measuring growth rates of deviations in carefully selected different directions. Typical 
deviations have a component in the most rapidly expanding direction, so they grow with 
the largest Lyapunov exponent. 

Of this chaotic behavior, little is seen on a macroscopic level, especially when we consider 
systems in or near equilibrium, but surprisingly, from numerical simulations, one found 
relations between Lyapunov exponents in stationary non-equilibrium systems, and linear 
transport coefficients Even if the generality of such relations can be questioned it 

is still interesting to consider models in which they hold. Among these are particle systems 
subject to a velocity independent external force and kept at a constant kinetic energy by 
means of a Gaussian thermostat 0J^. 

For the last six years, efforts have been made to get an analytic grasp on the Lyapunov 
exponents. For the Lorentz gas at low density in two and three dimensions. Van Beijeren, 
Dorfman and co-workers have set up a kinetic theory in which all Lyapunov exponents could 
be calculated and relations with transport coefficients could be verified p|-p!8| . The hard disk 



and the hard sphere gas, which were next in line for kinetic investigation, proved harder to 
deal with. It is possible to obtain estimates for the largest Lyapunov exponent in equilibrium 
for these systems, based on heuristic effective dynamics of deviations [p!^-pl|. Some results 



also have been obtained for the sum of all positive Lyapunov exponents, the KS-entropy 



|2T|,|22 



In these approaches, one uses the linear dependence of the deviations after a collision 
on their values before, together with the linear growth of position deviations during the 
times between collisions. From this, it was argued in Refs. ||l9| , |2l[| that at low densities, it is 
enough to look at the clock value 

\n{\Svi\/5vo) 

h = — — , (1) 

I lnn| 

where Svq is an arbitrary unit of infinitesimal velocity. These clock values approximately 
change in a collision between particle i and j according to 
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k[ = k'j = max{ki, kj) + l + 0{l/\nn), (2) 

which is vahd to leading order in n. The clock values k could therefore be restricted to integer 
values. In a chaotic system, each deviation dvi is expected to grow exponentially and the 
clock value ki to grow linearly with time. Because not all deviations of the particles are 
identical, there is a distribution of clock values. The dynamics described by Eq. (Q) tends to 
bring clock values far below average at any given time back towards the mean. As a result 
of this the distribution of all clock values with respect to the instantaneous average for 
long times approaches a stationary mean profile, about which the actual distribution keeps 
fluctuating. It then follows that the average clock value increases by 1 plus half the average 
difference between clock values per unit time (the mean free time between collisions for a 
single particle). In an abstract sense, this is equivalent to the situation of a propagating 
front, where a stable phase propagates into an unstable phase. The propagation occurs 
along One of the phases for a given "position" k is to have the clock distribution 

concentrated to its left, i.e., around lower clock values. This distribution will shift to the 
right and go beyond k. So this phase is unstable. The stable phase for a "position" k is 
to have the distribution concentrated to its right. On a technical note, it turns out that 



this problem falls into the class of pulled fronts [^, which is fortunate as it is known how 
to obtain the front propagation speed w for such systems. Using those techiques, for long 
times the average clock values was found to behave as 

N 

k = N'^ ^' W = ^0 + wut, 
1=1 

with u the average collision frequency and w a constant of order As a consequence of 
Eq. (|1]) the largest Lyapunov exponent is 

= —wulnn. (3) 

to leading order in the density. In this analysis was refined further. It was recognized 
that the distribution of clock values depends explicitly on the speed of the particles and the 
equations describing its time evolution were adjusted accordingly. Nonetheless the pulled 
front analogy remains fully valid. 

The purpose of this paper is to give a firmer basis to the heuristic arguments leading to 
equation (^, and to extend the methods so as to be applicable to higher densities and to 
non-equilibrium cases. 

A typical example of a non-equilibrium case is the hard disk gas in a uniform shearing 
state, described by the so-called SLLOD equations [^. The gas is confined between two 
infinite (one-dimensional) parallel plates a distance 2L apart, moving in opposite directions 
with velocities ±7^ (see Fig. jT]). We look at the limit of large L and infinite extension 
in the x-direction, with fixed shear rate 7 and density n = N/V. For small shear rate, a 
linear velocity profile develops in the system, so that the fiuid velocity at y is u = 'yyx. We 
remark that, although we want the volume V to be macroscopic, L should not be too large, 
otherwise the Reynolds number would become so high that the system becomes turbulent 
and the assumed linear velocity profile breaks down. 

The paper is set up as follows. In Sec. H we discuss the dynamics of deviations. In 
Sec. |Tl|we consider the velocity distribution function and the distribution of the duration of 
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intercollisional flights for tlie liard disk gas under sliear. In Sec. |^ we set up a generalized 
Boltzmann equation for a one particle distribution function that includes the deviations, 
and expand that in powers of 1/ In n. In Sec. |V|, this equation is reinterpreted as describing 
a propagating front. The largest Lyapunov exponent is proportional to the front velocity w, 
which we can determine using a standard method for pulled fronts ||2^. The perturbation 
expansion in the density is further developed in Sec. 0. In Sec. [VII| we calculate the first two 
terms in the density expansion of the largest Lyapunov exponent for the two-dimensional 
hard disk gas in equilibrium. The results are formally extended to the shear case in Sec. [VIII . 
We conclude with a discussion in Sec. [TX. 



II. CHAOS IN HARD DISK GASES 

The dynamics of the hard disk system is defined as follows. When disks i and j hit each 
other, they collide elastically. We define r^j = — Vj and Vij = Vi — Vj, and get 

V'. = Vi - {& ■ Vij)&, 

v'. = Vj + {& ■ Vij)&, (4) 

with <T the unit vector in the direction of the line connecting the center of the two disks at 
contact, i.e., a = a'^Vij. Primed quantities denote post-coUisional values throughout this 
paper. The positions remain unchanged, r'. = Vi, r'. = r 



In between collisions, the coordinates of disk / satisfy 

ri = vi ; vi = —Fi{{rh,Vh]), (5) 
m 

The forces Fi are smooth functions of the coordinates {r/i, d^}, h = 1 . . . N. In an instan- 
taneous collision between disks i and j, the smooth forces Fi and Fj cannot perform any 
action, therefore Eq. (^) describing the collision, holds for any Fi. 

The dynamics of the deviations 6ri and dvi in collisionless flight are given by the lin- 
earized version of Eq. (|) , 



dFi ^ OF, \ 

Srh,a + -Tr—5vh,a , (6) 




dVha 



where Sr^^a and Svh^a denote the a-th component of Svh and dv^, respectively. 

To find the collision dynamics, we use a method developed both by Gaspard and Dorf- 
man [EH] and by Dellago, Posch and Hoover ||26|. Here, we work out the dynamics for the 



general system of Eqs. (0) and (|^). For the equilibrium case this was done in Refs. |2T1 , |26 
The reference trajectory and the adjacent trajectory are infinitesimally close so we can as- 
sume that they have the same collision sequence. The most subtle ingredient in the derivation 
of the collision dynamics of deviations, is the time difference 5t between the collision 
on the two adjacent trajectories. We set the time of the collision equal to zero for the 
reference trajectory, so that of the adjacent trajectory equals 6t. We consider here the case 
that 6t is positive, but one easily checks that the final results are equally valid for negative 
St. We define 
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5ri = r*{0-) - n(O-), 6vi = v*{0-) - vi{0-), 
5r[ = r*{6t+) - ri{6t+), 5v[ = v*{6t+) - vi{6t+), 

where I = i or j, the superscript * denotes values on the adjacent trajectory and + and — 
indicate after and before colhsion, respectively. 

The time shift St can be found from the requirement that at the instant of collision, 
the two disks are a distance a apart, i.e., |T"ij(0)| = a and \r*j{St)\ = a. Because the time 
difference 6t is infinitesimal, we only have to express the r^{6t) to linear order in 6t, yielding 
r*{St) = T; + Sri + Vi6t. Note that here (and in the rest of this section) unprimed quantities 
without time specification are assumed to carry their value at t = before the collision, e.g., 
Vi = vi{0~). From the requirement |r^(5t)p — Ir^jp = 0, we get 



2rij ■ {drij + Vijdt) = 
to linear order in the deviations, so 



The difference in collision normal S& = &* — & follows from 

OC — — — / ^ s OVij, 

a a a[a ■ Vij) 

where we used Eq. (|^) in the last equality. 1 is the identity matrix and we use the conventions 
that non-dotted products of two vectors are dyadic products, and a product of matrices 
always implies matrix multiplication, as does the product of a matrix with a vector. 

Consider first the position deviations. For the reference trajectory we have, for I = i or 

J, 

ri{6t) = ri + v[6t, 

because the trajectory is determined by the velocity Vi{0~^) after the collision at t = 0. For 
the adjacent trajectory we have 

because this trajectory has velocity v^{0^), before colliding at time St. We write dr/ = 
r* + v^St — ri — v[St and insert the expressions for v[ from Eq. (^) and the one for St from 
Eq. (0), to find 

5r[ = dr, - {drij ■ &)&, 

Sr'. = drj + (Svij ■ &)& (8) 

(neglecting again expressions quadratic in deviations). For the velocity deviation vectors on 
the reference trajectory, one finds for I = i or j, 

vi{St+) = vi{0+) + -Fi ({rh,v'j) St, 
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Of course, only the velocities of the colliding disks, Vi and Vj, have really changed. We 

abbreviate Fi (^{vh, v^}) by and Fi {{vh, Vh}) by Fi from now on. For the adjacent 
trajectory. 



V 



m 



1 1' 
v? + —Fi5t 

m 



where in the last equation, the difference between F^5t and Fi5t has been neglected (second 
order in the deviations). The prime means that the collision rule for velocities should be 
applied. For the adjacent trajectory, this is a collision with collision normal ct*, so, for disk i, 



v*{6t-)-a* [&*-v*^{6t- 
Vi{6t^) + dvi — & {& ■ Svij) 
-[(6- ■ Vij)l + &Vij]d& 

+-{F, - Fl)6t - 1 . {F, - F,}) a6t. 

m ^ m 



Inserting the expressions for 5& and 5t, we obtain 



6v[ = 5vi — a-{a- ■ Svij) — (Q^ — £i)Sr, 



1-3 



where 



Si 



{{a- ■ Vij)l + &Vij\[{a- ■ Vij)l- VijO-] 

a{& ■ Vij) 
[a.{F,-F,)]a + F:~F, a 



(9) 



m 



(T ■ Vi 



Similarly, for disk j, we find 



dvj + &(& ■ Svij) + (Qa- - Si)dri 



Thus when the dynamics of a point in phase space is determined by Eqs. 
dynamics of the deviation vectors is given by Eqs. (||) and (p|-p!0|). 



and 



(10) 

, the 



III. SHEARED SYSTEM 



In the sheared hard sphere gas, it is convenient to transform to the peculiar velocities of 
the particles, 

Vi = vi- u{ri) = vi- -fyix. (11) 
The equations of motion in coUisionless flight (r; = Vi, vi = 0) are transformed to 



ri = Vi + jyix, 



(12) 
(13) 
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These equations are called the SLLOD equations of motion [^]. The pseudo-force —m'yVi^yX 
is called the shear force. Due to Eq. (|l3), in a collision between disks i and j, the peculiar 
velocities transform as 

Vf = Vi- {a- - Vij)d- - a-fd-{x ■ &){y ■ &), 

V; = V, + (a ■ V,j)a- + a^&ix ■a)iy- a), (14) 

In simulations, boundary effects can be minimized by a special kind of periodic boundary 
conditions, Lees-Edwards boundary conditions, in which case the periodic copies in the y- 
direction are moving with respect to another. When a particle crosses the boundary, it is 
put back at the other end but with a corrected position and unchanged peculiar velocity 

iJUJH- 

The system as defined above does not have a steady state. The reason is that the 
shear force continuously converts macroscopic kinetic energy of flow into heat, i.e. internal 
kinetic energy. In realistic Couette-flow situations the work required for this is performed 
by the shearing walls of the system. In the present situation this work is a consequence 
of the boundary conditions. Because the location of the y-coordinate of the Lees-Edwards 
boundary is arbitrary, the system develops no temperature gradient, in contrast to a system 
with realistic boundaries. In such a realistic system a stationary state is usually reached by 
the establishment of a stationary heat flow towards the boundaries, which absorb the heat 
and transmit it to a thermostat. If one doesn't want to include the environment explicitly, 
as is usually the case when one performs MD-simulations, one needs to put in a mechanism 
by hand to extract heat from the system. Such a mechanism commonly is also called a 
thermostat. Several thermostats are around for non-equilibrium systems, but we focus on 
one in particular. An extra term is added to the equations of motion for the velocities during 
coUisionless flight, which become 

V, = --iVy^x - aVi. (15) 

We choose for a a constant positive value, so that the extra terms can be interpreted as 
standard friction forces. In molecular dynamics simulations it is more common to use an 
isokinetic Gaussian thermostat, which keeps the (peculiar) kinetic energy 



strictly constant. In that case, a depends on the positions and velocities of all the particles 
and can be chosen such that the equations of motion are time-reversible, i.e., form invariant 
under a time reversal operation, a then may take both positive and negative values and 
only on average will it be positive in the stationary state. In the thermodynamical limit, 
this thermostat is equivalent to the one with a constant a, p8| , |29[] and we choose the latter 
one, as it is simpler. The value of a determines the average peculiar kinetic energy, which 
is identifled with the steady state temperature through 

The brackets here denote a time average, which, when the system is ergodic, is also the 
average over an appropriate steady state distribution. Our thermostat suppresses turbulence 
p3| in regimes where it is expected physically. So the model is representative for an actual 



physical system only for low enough Reynolds number. 
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In the low density regime we can use the Boltzmann equation for the one-particle distri- 
bution function |1^.3^. We only consider stationary flows of uniform shear, so we can focus 



on the velocity distribution function of the form f{V), which is normalized to unity. The 
Boltzmann equation for this distribution is 

-dv^-[i7V,yX + aVi)h] 

= j ■■■ J'na\&-Vu\{M-fif2)dV2d&, (17) 

with the short hand notation /j = /(V^) and f- = f{Vi'). The prime on the integration 
denotes the condition & ■ Vu < 0. 

We want to discuss briefly how one can solve this Boltzmann equation for small shear 
rates. In equilibrium, 7 = a = 0, and the right-hand side vanishes if / is a Maxwell 
velocity distribution. To allow for treating the left-hand side as a perturbation, 7 has to be 
small compared to ?T,a(|i;i2|). Hence, the small parameter proportional to the shear rate is 
7 = 7a/(nfo). A Chapman-Enskog expansion of Eq. (^) can now be made by expanding 
in powers of this parameter. The expansion of / has the form 

f{V) = ^{V)[l + + fg^iV) + ...], (18) 

where (p is the Maxwell distribution 

In the rest of the paper, we use 7 as an independent variable, separate from the density 
h. In this spirit, higher density corrections are terms which vanish (relatively) when the 
density is lowered while 7 is kept fixed. Thus we neglect the last terms ±07(0* ■ x){& ■ y) in 
Eq. (p!4D, as they are (9(7nfo), i.e., one order in density higher than the other terms, which 
are of O{vo). Likewise, the restriction on the integral in Eq. (|17D is replaced by Vu ■ o- < 
and u{a&) is neglected in |(t ■ Vi2\ = \a- ■ {V12 + u{a&))\. Finally, Vij is replaced by Vy in 
the matrix Q^. in Eq. (|^). 

Within the framework of the Chapman-Enskog expansion the zeroth order solution, (p, 
entirely determines the temperature T; higher orders should not change the second moment 
of /. For fixed 7 and T the friction coefficient a has to be expanded as a = 7[ai + 702 + •••]• 
To determine the coefficients, we note that {dtJ2i\Vi\'^) = in the stationary state and that 
X]j is conserved in a collision, so that for the time derivative we can insert Eq. (plsl). In 
this way, we see that 

« = 7{V.Vy)/{\V\') = ^J V.VyfiV) dV. (20) 

Thus, if we know / up to order 7""^, we can calculate a up to order 7", which we need to 
calculate / up to order 7". One also immediately sees that ai = 0, because the average of 
the odd function VxVy with the even distribution ip is zero. This is no surprise, as a gives 
the dissipation and this should be an even function of 7. 

We will need the distribution s{V, r; t) of particles at time t with velocity V and time 
T passed since their last collision. In the stationary state this satisfies the equation 
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dy ■ [V s( V, r)] + drs{V, r) = -v{V)s{V, r) 
for r > 0. The collision frequency v is given by 

v{Vi) = j ■■■ J'na\a ■ Vi2\f2dV2d&. 



(21) 



We rewrite the density s in terms of a conditional distribution function: s{V,tau) = 
f{V)S{T\V). With 

u*{V) = u{V) + dyV + V-dY\nf, (22) 

we find the following equation for S{t\V), 

drSirlV) + V ■ dySirlV) = -z/*(V)5(r| V), 

with the initial condition to be determined by normalization. The general solution of this 
equation is 



5(r|y) = So(y(-r))exp 



i^*iVit'))dt' 



where V{t) is the solution of the equations of motion, Eq. ([T5|) , with initial condition V{0) = 
V, i.e., V(t) = e~°*[V — xjtVy]. It is simple to show that for /q°° S{T\V)dT to be equal to 
1, we need So{V) = z/*(y), so 



S{t\V) = z/*(V(-r))exp 



u*{Vit'))dt' 



(23) 



IV. GENERALIZED BOLTZMANN EQUATION 

In this section we derive a generalized Boltzmann equation for the distribution function 
of V and 6V. But first of all, we make the dynamics of the deviations explicit for the case 
of the SLLOD equations (^^-|TB|) that we are investigating. In terms of rj and Vi, we have 
Vi = a'yyiX — aVi. For this case, we get 

a ■ X a ■ y ^ ^ 

Ci = a'ya crcr, 

& ■ Vij 

(see Eq. The SV then change in a collision according to 

SV^ = SV,- &{& ■ SV,,) - (Q^ -S,- S2)5n, 

6V; = 5Vj + &{& ■ SVij) + (Q^ -Si- S2)6rij (24) 

where S2 = 'y[x{& ■ y)& — a-{x ■ o-)y]. The terms Si and S2 are of higher order in the density 
than Qo-, which is of order (fo/a), as one sees when expressing these quantities in terms 
of 7 and h; Si = O^VQ^^n^ /a) and S2 = O{vo'jn/a). As we are restricting ourselves to the 
leading powers in density, with at most correction terms of relative order 1/| ln?7,|, we neglect 
Si and S2 in the sequel. 
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In collisionless flight, the dynamics of deviations is, 



The solutions of these equations are: 

5V,{t) = e-"*[5F,(0) - x^tdVyM] (25) 

1 - e""* 

5ri{t) = 5r,,(0) + 5yi{Q)-ftx + 8Vi{Q) 

a 

+ — -o l5Vi^y{^)x 

= tStSV.iO) + Sr.iO) + 5yatx. (26) 

We follow the deviation dynamics from collision to collision. So we are interested in the 
case where the time t in the equations above is of the order of a mean intercoUisional flight 
time, t = 0{a/ iyofi)). Then at = 0(7^), and •yt = 0(7), and both are of zeroth order in 
the density. We assume that 

OiSVi) = 0{dn) X vo/a, (27) 



just after a collision. Then, at the next collision, the terms with a factor tSVi{0) in Eq. (p6D 
typically are one order in 1/fi larger than the corresponding terms proportional to <5rj(0), 
which therefore may be neglected. The Svi deviations just before a collision then are 
<5Tj(rj) = TiSriSVi{0), where Tj is the time from the previous collision of particle i. The 
effective equation linking the dVi just after collision to their values just after the previous 
collisions of the two particles, is found from Eqs. (^4]) and (^): 

6V^ ^ -(5V; ^ -Q^{r,Sr,SV, - r.Sr^dV,), (28) 

where we neglected terms of higher order in n coming from the center-of-mass contribution. 
To check for the consistency of our assumption about the relative order of velocity and 
position deviations just after a collision, we consider also Svi and Svj after the next collision: 

6r[ ^ -dr'^ ^ (il - &&){rSnSV, - TjSr^SVj), 

and we see that if Eq. (|27|) holds for Svi, dVi and dvj, 6Vj just after the previous collision, 
it also holds for Sr' 6V/ and 6r' SV-. 

In two dimensions, the matrix Q^- in Eq. ( pSD is 



a(<T ■ Vij) 



where R denotes a rotation over 90 degrees counterclockwise. To leading order in the density 
we may replace by Vij. Neglecting further Si and S2 in Eq. we obtain as effective 
equations of change for the velocity deviation vectors in a collision. 



SV[ = \i \^ . ^, (29) 



a{& ■ Vij J 



dV'^ = -6V[. 
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To compare different contributions to dV[ and SV'j, we want to know the order of dVi and 
SVj. We write the value of dVi just after the previous collision as 



6Vi = vo 



n 



VqC 



(30) 



where e,- is a unit vector and is defined as 



Innl 



The clock values and unit vectors thus defined do not change during collisionless flight, only 
in collisions. In contrast to earlier work in Refs. [0,^ , the clock values ki are real numbers 
here, not integers. 

We obtain k[ from Eqs. (EPI) and 



k[ = ^In 



a\a-v 



21 



and distinguish two cases, telling us which of the two terms inside the logarithm is the 
largest, and hence, the most important. We define 



b{V, r, e, &) = In 



rnKV ■ S^e 



aV ■& 



(31) 



(32) 



and say that 

/: disk 1 dominates disk 2 if 

ki + ^96(Vi2, Ti, ei, &)> k2 + i^> b{Vi2, T2, es, 6") 
and that 

//: disk 2 dominates disk 1 otherwise. 

We also define an "alignment" criterion: 

(+) The unit vectors ei and 62 are said to be aligned if RVl2 ■ 8^161 has the same sign as 
RV12 ■ 8^-262, 

(— ) and anti- aligned if they have opposite signs. 

From now on, when a ± or =p occurs in an equation, the upper sign corresponds to the 
aligned case, the lower one to the anti-aligned case. With these definitions, we write for the 
case /: 



k[ = ki + l + i^b{Vi2,n,ei,&) 

jk,-k^)/i) ^MKVi2,T2,e2,a)] 



+t91n 



l^e^ 



exp[6(Vi2,ri,ei,6-)] 



(33) 



and for case 11: 
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k[ = k2 + l + ^96(Vi2, T2, 62, a) 



(fc,-fc,)/^ exp[&(Vi2,ri,ei,o-)] 
exp[6(Vi2,r2,e2,o-)] 



(34) 



At this point it can be made clear why the distinction of one disk dominating over another 
was made. Consider Eq. (^3|). Because is small for low densities, and because Eq. (|32D 
holds, the term after "l=p" inside the logarithm tends to be small, at least if ki and k2 
differ by an amount of 0(1). Consequently, this whole term is small, or, more precisely, it is 
appreciable only for \k2 — ki \ = 0{'d). Also the term -db is typically small. Hence k[ = ki + 1 
almost always, in the limit that (the infinitely dilute gas). The same limit in case // 
yields k[ = k2 + l. We see that indeed, the "dominant" particle determines the value of the 
clocks after collision, at least if the density is low enough. This limiting dynamics for low 
density, expressed in Eq. (g), was derived before in Refs. [^] and ||2T| and proved sufficient 
for obtaining the leading term in the density expansion of the largest Lyapunov exponent. 

We consider the conditional probability distribution function for having clock value k 
and unit vector e just after a collision at time t, given that the post-coUisional velocity is 
V. This function is denoted by /(/c, e\V; t) and it obeys 

u:{k,e,V-t)=u{y)f{k,e\V-t), (35) 

where uj{X\t) stands for the rate at which particles with attributes X are produced in 
collisions at time t. The production rate u}{V) of V in the stationary state is independent 
of time and satisfies 



Vi) = j J'an\& ■ Vu\f[f^dV2d& = z/*(Vi)/i, (36) 



where we used Eqs. (|15|) , (p!7|), (^11) and (|22|) . By considering the rate of restituting collisions 
that produce the right {k, e, V; t), we find that uj{k, e, V; t) satisfies the equation 

u{k,e,V-t) = j ■■■ f an\& ■V^2\flf2SiT,\V,)Sir2\V2) 

xf{ki,ei\Vi{-Ti);t-Ti) 
xf{k2,e2\V2{-T2);t- rs) 
x6{V[ - V)5{k[ - k)6{e[ - e) 

xdVidV2dkidk2deide2dTidT2d&. (37) 

In the arguments of / the velocities need to be traced back to the previous collision, because 
/ was defined in terms of the variables at that instant of time. Hence the appearance 
of V{—t). The clock values and unit vectors do not need such a correction, as they do 
not change in between collisions. We symmetrize the equation with respect to e, because 
e — *• — e only means interchanging the reference and the adjacent trajectory, and this cannot 
affect their rate of separation. Hence we can replace S{e[ — e) = 6(RVi2 ~ ^) by [5(RVi2 ~ 
e) + 5{RV[2 + e)]/2 = ^S{e-v[,). 

Through linear order in the logarithmic terms in the expressions (P3|) and (53) for k'l 



may be ignored in Eq. (|3^); their inclusion merely gives rise to corrections of 0{i)'^). With 
this approximation Eq. (p^) may be rewritten as 
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u!{k, e, V; t) 

= ^/--- f S{n\V,)S{T2\V2)fj2 

xan\a ■ Vi2\6{V[ - V)i(5(e ■ v[^) 

X C( A; - 1 - Vi2 , ri , ei , 6-) , ei I Vi (-n ) ; t - n) 

xC{k-l- m{Vu, T2, 62, 6-), eal V2(-r2); t - T2) 

xdVidV2deide2dTidT2da- + 0{^^). (38) 



Here, we introduced a cumulative distribution, defined as 

rk 



C{k,e\V;t)= f f{k*,e\V;t)dk*. (39) 

J— CO 



Integrating Eq. (^) from —00 to k, using Eqs. (^), (P^ ) and (PD| ) for the left hand side, 
and changing integration variables from precoUisional velocities to post-coUisional ones, we 
arrive at 

ly*iV^)C{k,e\V,■t) 

= J ■■■J's{n\V[)S{r2\V',)fr'f[f!, 

xC(A; - 1 - m{V'^,, n, ei, &), ei|y;(-ri); t - n) 

xCik-1- ^biV^',, T2, 62, a), 62! V'2(-r2); t - T2) 

xan\a- ■ Vi2||5(e ■ V i2)d'V2deide.2dTidT2da- . (40) 

This is the generalized Boltzmann equation for the cumulative distribution function C, up 
to 0{d'^), which will be the starting point for our calculations of the maximal Lyapunov 
exponent. 



V. FRONT PROPAGATION 



The generalized Boltzmann equation (^D|) can be interpreted as describing a propagating 
front. The propagation here occurs on the real line of clock values k: As clock values tend to 
grow there is a movement towards higher clock values. The two "phases" that are separated 
by the front are the stationary solutions C{k,e\V]t) = on the left (no particles have a 
clock value smaller than the fc-values in this region) and C(fc, e|; V;t) = P{e.\V) on the 
right (all particles have a clock value smaller than the /c- values in this region). P{e\V) is 
the conditional probability for a particle to have unit vector e given that its velocity is V, 
i.e. 

P(e|Vi)= / f ^^^na\&-Vi2\\K^-Vi2)dV2da. 

It is easy to see that C = is stable, whereas C = P{e\V) is unstable. In the simplest 
situation, the front has a fixed shape and moves to the right with a constant velocity: 

C{k,e\V,t) = F{x,e\V) (41) 
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where x = k — wut with u the average colhsion frequency as mentioned already in the 
introduction. The constant w is called the clock speed. For the simpler clock model where 
each particle is fully characterized by a single clock variable which increases with time 
according to Eq. (H), it turned out that the front falls in the class of so-called pulled 



fronts, as opposed to pushed fronts |2^. We will assume that in the present model this is 
also the case. 

In short, the asymptotic front speed of pulled fronts is determined as follows. Insert a 
propagating front solution like Eq. (|41|) into the front equation (|40|), and linearize around 
the unstable phase. The resulting linear equation has solutions which are linear superposi- 
tions of exponential functions of the form e~^^ (multiplied by a polynomial in x in case of 
degeneracies). For fixed w, there are a number of possible values of s (typically infinite but 
countable). The dominant term in the superposition is the one where s = Sd has the small- 
est real part. Since C has to be monotonic in x, the asymptotic large x behavior ~ e"^"^^ 
has to be so too, hence Sd has to be real. This turns out to be possible only for w larger 
than some critical value w*. So the asymptotic speed, if it exists, is greater than or equal 
to w*. In fact, for a large class of initial conditions, the asymptotic clock speed is exactly 
w*. Especially, this is true for localized initial conditions, implying that all initial fc- values 
fall within a finite range (or in fact, are smaller than some value kmax] localization on the 
small-fc side is not important). The same speed is also selected by initial distributions that 
fall off sufficiently rapidly (typically faster than any exponential of type exp[— ca;]) at the 
large-fc side. For systems with finite numbers of particles, which we are typically interested 
in, the initial distribution is always localized, leading to automatic selection of the minimal 
clock speed w*. The additional effects of finite particle numbers only reduce the clock speed 
further, because correlations between particles tend to reduce their clock speed differences 
and thereby the boost a "slow" particle receives from colliding with a "fast" one A 
more detailed exposition of the velocity selection and other aspects of pulled and pushed 
fronts can be found in Ref. | |24|| . 

Applying this scheme to the generalized Boltzmann equation Eq. (0), we first insert the 
Ansatz (^Tf) and consider the resulting equation. We find 

u*{V^)F{x,e\V,) 

' S{n\V[)S{T2\V',)f^'fJi 

xF{x-l-^b{Vi2',n,ei,&)+wun,ei\V[{-Ti)) 
xF{x-l- i^biVu', r2, 62, 6") + wur2, V2(-^2)) 
xan\a- ■ Vi2\^6{e ■ Vi2)dV2deide2dTidT2da-. 

We linearize this equation writing F{x, e\V) = P(e| Vi) — A(x, e| V). The resulting linear 
equation for A has superpositions of exponentials as solutions. It turns out convenient to 
represent these as 

A(x,e|V) = E^^^^^^^^^Me,V)e--^/-. (42) 

The characteristic values Si and corresponding characteristic functions Ai are solutions of 
the linearized equation with A taking the form 
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This gives the following characteristic equation 
K{v*{Vr) + sv)A{e,Vr) 



X Ps(r|K)exp[-(lnA)^6(V'2i,r,e*,^) 

k=l,2 

* -17-1 



X A(e*, Vk{-T))dV2dedTd&, 

where we defined the eigenvalue A = e~^/^ and 

P^{r\V) = {v*{V{-T)) + sv) 

rO 



xexpl-y {u*{V{t)) + su)dt 

We make one more expansion in -(9: 

A{u* + su)A = L°i - t^lnAL^i, 
where we left out the argument of u*. The operators L° and are defined by 



f[f^f,-'na\a■V^2\^5{e■V^2) 



X J2 Ps{T\V',,)A{e\Vi{-T))dT*de*dV2dTda, 

k=l,2 

1 



X J2 Ps{r\V',,)b{V2i,T,e*,&)A{e,V'j,{-T))dT*dV2dedTda. 



f[f^f,-'na\&■V^2\-6ie■V^2) 



fc=l,2 



(43) 



(44) 



(45) 



(46) 



Finally, Eq. (0) can be transformed into an equation in which only the integrated function 
over e enters, i.e., A{V) = [PA]{V) = {2n)-^ J de A{e,V). Defining 



{e,Vi) = J---J f[f2fina\a ■ Vi2|vr5(e • V12) 
X Ps{T\V'^)A{V'k{-T))dV2dTda, 

k=l,2 



we see from Eq. that we can write 



_ f 



(47) 



We see from Eqs. ( ^ and ( ^ that A = A ^{u* + su) ^L°y4 up to first order in {}, so 
Eq. ( PI) can be written as 

A(z/* + siy)A = L°A - ^9i^L(z/* + sz/)~^L°A + 0{^^). 



15 



Applying P yields a closed equation for A: 

A{u* + sij)A = L^A - M"^ In AL^A, (48) 

where 

L° = PL°, 

and 

= PLi(z/* + sz/)-^L°. (49) 

These are just the first few operators that appear in an expansion in t?, each still containing 
all orders of the shear rate 7. We use Eq. ( ^ ) to find A as a function of s through linear 
order in ■(9. As A will be equated to e~^/'^, we are interested in the largest eigenvalue A(s) 
for which Eq. ( |I8| ) can be satisfied, because this corresponds to the most slowly decaying 
mode in Eq. (0). There are real solutions to A(s) = e"'^^^ for s if w > w*, for w < w* 
only complex solutions exist, leading to undesirable oscillations. To find w* we look at the 
function 

There are no real solutions of w{s) = w if w is smaller than the minimum of the function 
w, hence this minimum is w*. 



VI. PERTURBATION IN THE DENSITY 

This section will be devoted to finding the solution of Eq. (^) in a perturbation expansion 
in that is we will find A(s) and with that determine the minimum of — s/lnA(s), which 
is the clock speed w. 

We consider Eq. (|48|) for the largest eigenvalue A(s). Let us assume that Eq. (^8[) is 
solved to zeroth order by A° and A°(s), i.e. 

A°(s)(z/* + sz/)A° = LM°. (51) 

A^ depends on s, although we will not denote this explicitly. We remark that this zeroth 
order equation was solved in Ref. [^] for the equilibrium case. In the general case, L° is 
not self-adjoint, so we need the left eigenfunction too. We denote this function by A^. An 
inner product is defined as 

{A,B) = J A{V)B{VMV)dV, (52) 

where the Maxwell distribution was given in Eq. (p^). The function A^ may be chosen such 
that under this inner product {A^, A^) = 1. Inserting the expansions A = A^ + i)A^ + 0(i?^) 
and A(s) = A°(s) + '&A^{s) + 0{'d'^) into Eq. (H3) and considering the 0{'d) terms, we get 



16 



A"(s)(z/* + sp)A' + A'(s)(z/* + sp)A 

The inner product of this equation with yields 

^ lnA°(s) (A°,LM°) 



(53) 



where we used Eq. The critical value w is found from {dw / ds){s*) = and w = w{s*) 
(where s* is the location of the minimum). Using Eq. (|53D for A^(s) and Eq. (|50| ) we find 



Wis) = w 



7,0 1 



s) + {}w'^{s) + 0(^9^), where 



InAO(s) 

7,0/ 



w°(s) (AQ,LM°) 
AO(s)M^°,(i^* + si^)^°)' 



(54) 



Note that A° and depend on s as well. Again we assume the minimum of to be 
known, and to be located at s°, with a value of = w^{s^). The first order value of the 
minimum, which to first order is located at s* = s° + "^s^, is (up to that order) 



7,0 



dw' 
1^ 



where we used that the derivative of at s° is zero. From Eq. (|5^ and the identity 
A" = exp{—s^/w^), the value of the critical clock speed up to order follows as 



w = w 



1 + ^e^^V-" (v4°,LiA°) 2 

^ {AO,iu* + s<'u)A^)^ ^ 



This is the value of the clock speed that enters into the Lyapunov exponent, according to 
Eq. (D, so 



In — + e 

n 



(AO,LiAO) 



(AO, {v* + s^u)A^) 



(55) 



Correction terms to this expression are 0{vd). 



VII. EQUILIBRIUM CASE 

We will now explicitly evaluate Eq. (^) for the equilibrium case. In that case, i.e., 7 = 0, 
Eqs. (^) and (^) simplify strongly. One simplification is that the velocity distribution 
is known to be the Maxwellian, i.e., f{V) = v^(V), so in the integrand, /i'Vi/2 = '^2- 
Furthermore V = 0, which, according to Eq. (P^D, makes i'*{V) and i^(V) equal to the 
equilibrium collision frequency i^oiV) of a particle with velocity V. The expression in 
Eq. (^3D now gives 
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Ps{t\V) = iMV) + su) exp {-h(V) + suo]t} . 
The expression for the function b in Eq. simphfies because St = 1: 



Tfl 


RV 


e 


a 


V a 





b{V, r, e, <t) = In 



We see that now the integration over r in Eqs. ( ^31) and (|^) can be performed. The result 
is that Eq. (^81) becomes 



^^0 



(56) 



where 



(Vi) = J ■■■J na\& ■ Vi2\ [A[ + A'^] ^2^^26/6", 
Here A'^ = A{V'i^), and, analogously to Eq. (|49|) , 



where 



na\& ■ V12I ^ ttA'^ 

k=l,2 



Ua 



x6{e ■ Vi2)(f2dV2d& 
r' 1 



-na 



X In 



2 



V12 



IRV12 ■ e* 



MVk) + si^o \Vi2-o-\ 
x6{e ■ V i2)ip2dV2d&de* 

Here C is Euler's number, 0.577. . .. In general, subscripts denote the equilibrium values 
of quantities introduced before. 

To zeroth order in i), Eq. (^) reads 

A°(z/o + sz>o)AO = L'X- 

This is the very same eigenvalue problem that was found in a more heuristic derivation of a 
Boltzmann equation for clock values in Refs. [^, with exactly the same operators, except 
that s was represented as ■jw and Lq and 1^0 + suq were denoted as L and i^oWg, respectively. 
The largest eigenvalue Aq(s) was determined numerically in that paper, for a range of s, such 
that the minimum of w could be determined. Briefly, the method used was the following. 
A basis of functions was constructed that are orthogonal with respect to the inner product 
defined in Eq. (H), starting with 1, V /vq, l\V/vo\'^ - 1 and l\V/vo\'^ - |V/woP + 1- Finite 
matrices were constructed containing the matrix elements of the operators with respect to 
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the four mentioned basis-functions. These finite matrices were used in Eq. (|56|) instead of 
the real operators, to get a numerically feasible eigenvalue problem. In Ref. |2l[] it was 
checked that omitting the basis function containing only changes the result by a few 

tenth of percents, so that is the accuracy of the numerical results. In this paper we choose 
to work with just 1, V/fo, and ||V/fop — 1, as this simplifies the calculations. On that 
truncated basis the value and position of the minimum are found to be 

w° ^ 4.732, (57) 

s° ^ 3.506. (58) 

and the eigenfunction Aq is 

.4°(y) ^ 0.612 + 0.194i;(72|v|2_ (59) 

The value of Wq gives the leading behavior of the largest Lyapunov exponent: Xq = 
-w{Ji^olnn + 0(l). 

Using these results we can obtain the first correction term. For that, we adapt Eq. ( |55D 
to this case: 



A+ = WqUq 



Note that Lq is self adjoint and therefore the left and right eigenvectors are the same. 
We have calculated the matrix elements {Aq, LqAq) and {Aq, {uq + SQZ/o)v4g) numerically 
respectively analytically, using the numbers of Eqs. (0), (|59|) and {^7\). The numerical 
integration uses a sampling from the Maxwellian and subsequent averaging of the rest of the 
integrand appearing in the matrix elements (including the collision frequency, for which we 
used a numerical approximation). The results are (AQ,LQy4Q) = — 10.85 nawo and {Aq, {uq + 
SqI/o)^o) = 19.28 nafo Thus we obtain the result that 

A([ = 4.732z/o[- Inn -2.48 + 0(1/ Inn)]. (60) 



VIII. PERTURBATION IN THE SHEAR RATE 

The result in the previous section can be the starting point of a perturbation theory 
in the shear rate in the case that 7 is non-zero, but small. From symmetry (7 — > —7), it 
follows that the Lyapunov exponent A"*" is even in 7, so to see the effect of the shear we 
need a second order perturbation theory at least. We will sketch the solution formally in 
this section, and leave the explicit calculations for the future. 

We start with the eigenvalue equation in Eq. ( ^8] ) to zeroth order, i.e., 

L°A° = A°(z/* + sz/)A°, 

and write = L[] + 7L0 + ^^Ll + 0{^^), u* = uq + ^vl + I'^v^ + 0(7^), = + 7^? + 
72^0 + 0(7^) and 

A%s) = Al + ^'Al{s) + 0{f) (61) 
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where we used that A'' is an even function of 7. For the same reason we may expand v as 
^ = + Substituting these expressions into the eigenvalue problem and equating equal 
powers of 7 one finds 



4° 

^05 



A\ = n [A°z/* - L? 
_ {Al [L° - A°(^* + s^,)]Al + (L° - A.lul)A\) 
{Al (z/o + suo)Al) 
Al = 7^(A°z/* - L?)A? 

+7^[A0(z/* + ^2s) + AO(z/o + sz7o) - Vi]Al 

where TZ is the inverse of Lq — Aq(z/o + sz/q), restricted to the subspace orthogonal to Aq (so 
A\ and A^ are made unique by requiring orthogonality to A^). 

From Eqs. (pO|) and (|61|), we obtain ^"(5) = Wo(s) + 7^^2(3) + 0(7"^), where 

^°(s) = -s/lnA°(.) 

The location of the minimum of the function w^{s) is shifted by an amount of order 7^, in 
fact 

dw'i, / ON 

^° = ^S-fli^ + 0(7') (62) 

as can be found from {d/ ds)w^{s^) = and {d/ds)wQ{sQ) = 0. Because we are expanding 
around a minimum, this shift is not needed for the value of the minimum of the function 
w°(s), i.e., 

Thus we can determine the shear corrections to the leading density term in Eq. (^). 

To calculate the density correction in Eq. ( |^ we do need the shift of the minimum as 
it enters in the matrix elements. This means that one has to insert Eq. (|62D and expand 
all relevant matrix elements in powers of 7. All of these can be obtained with the help of 
Eqs. (|18D, (13), (IID, (H), (H), (H), (H), (|31|), (H) and (||). In the end the evaluation 
will have to be done numerically. 



IX. DISCUSSION 

In this paper we developed an analytic method for calculating the largest Lyapunov ex- 
ponent of a many body system, based on the microscopic equations of motion. To be specific 
we restricted ourselves to uniform hard disk systems in two dimensions at low densities, but 
not necessarily in equilibrium. As a particular case we considered the uniformly sheared 
hard disk gas. 

To obtain the largest Lyapunov exponent we derived a generalized Boltzmann equation 
that describes the time evolution of a distribution function of particle positions and velocities. 
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together with deviation vectors in tangent space. At low densities the position deviations 
turn out to be unimportant, to leading orders in the density. The velocity deviation of a 
particle may be represented conveniently in terms of the logarithms of its norm- the so-called 
clock value k (Eq. (|I])) - and an additional unit vector, which in the end plays no essential 
role. 

As in the preceding papers |[T9|-pT| the generalized Boltzmann equation may be reinter- 
preted as describing the propagation of a pulled front on the real line of possible clock values. 
Therefore, by standard techniques the linearized version of this equation may be used to 
obtain the asymptotic speed of propagation of the front, which is directly proportional to 
the largest Lyapunov exponent. 

A remarkable property of the generalized Boltzmann equation is that its natural density 
expansion does not proceed in powers of the dimensionless density n, but rather in powers 



of ^9 = 1/1 lnn|. To lowest order, this reproduces the results of earlier work [jlOl , ^ , which 



therefore finds a firmer basis here. Here we also give values for the next order in 'd con- 
tribution to the largest Lyapunov exponent, and we give explicit expressions for the first 
non-vanishing shear rate dependent contribution to this exponent (quadratic in the shear 
rate) in the uniformly sheared system. 

In we made comparisons between our results and the results of numerical simulations 
of hard disk Lyapunov exponents by Dellago and Posch [2^. Within the numerical accuracy 
good agreement was found, but it turned out there are two complicating factors. The first 
one is that the density should be very low to be in the regime where a power series in 
1 / ln(n) can be expected to converge rapidly. The second one is that there are large finite 
size effects for the front speed w. One cannot (yet) reach the required number of particles 
in simulations for these effects to become negligible. The finite size effects can be estimated 
for asymptotically large using front propagation techniques, and they would scale as 
l/ln'^{N). It is not known at which number of particles this asymptotic result suffices, 
though. 

To conclude we mention some possible extensions of this work and some interesting prob- 
lems that are still open for research. First of all the expressions for the shear rate dependent 
contribution should be worked out numerically. We plan to do this on short notice. Then one 
would like to go more general potentials, to general densities and to three-dimensional sys- 
tems as well as two-dimensional ones. We expect that, as long as one stays at low densities, 
the generalization to other simple (short ranged and spherically symmetic) potentials should 
not pose any serious problems; one just has to use the generalized Boltzmann equation that 
is appropriate for the potential under consideration. Similarly we expect the generalization 
to three dimensions just will cause some additional technical complications, related to the 
fact that the directions of the velocity deviations after a collision will depend on those before, 
unlike in two dimensions. Generalizations to higher density will be much harder to accom- 
plish. Already the fact that the natural expansion in density proceeds in powers of i) rather 
than h indicates that a systematic expansion up to appreciable density will be forbiddingly 
hard. On the other hand non-systematic approaches, such as a generalized Enskog equation 
for hard disks or spheres, may give good approximations, but this has not been explored 
yet. 

One can also try to extend the methods developed here so as to calculate additional 
Lyapunov exponents. For non-equilibrium systems the most negative Lyapunov exponent is 
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For the SLLOD 
there are deviations of this rule 



of primary interest, because the sum of this and the largest Lyapunov exponent is directly 
related to the shear viscosity through the conjugate pairing rule ||2^ . 
equations with a Gaussian or "constant-a" thermostat, 
|3T|^3|, but these are of higher order in the shear rate, so that the connection with the 
linear viscosity still stands. 

Finally, at low densities, the methods developed here may probably be combined with 
those of Refs. []T2|,p2[ in order to calculate the sum of all positive Lyapunov exponents, or 
the KS-entropy, for non-equilibrium cases. However, a theoretical calculation of the full 
spectrum of Lyapunov exponents, which probably requires similar techniques, remains a 
very challenging open problem. 
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FIGURES 

FIG. 1. Velocity profile in a gas under shear 
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